/*
 * Equilibrator.cpp
 *
 *  Created on: 2 Aug 2011
 *      Author: Allan
 */

#include "Equilibrator.h"

// C++ includes
#include <algorithm>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <set>
#include <tuple>

// GeoReact includes
#include "Multiphase.h"
#include "ReactionSystem.h"
#include "Types.h"
#include "Utilities.h"

#include "Utilities/Algorithms.h"
#include "Utilities/PhaseUtils.h"
#include "Utilities/ReactionUtils.h"
#include "Utilities/SpeciesUtils.h"
#include "Utilities/StringUtils.h"

Equilibrator::Equilibrator()
{}

Equilibrator::Equilibrator(const System& system) :
system(system), eSolver(system.GetMultiphase(), system.GetEquilibriumReactions())
{
	// The number of species in the system
	const unsigned N = system.GetMultiphase().GetSpecies().size();
	
	n = VectorXd::Constant(N, 1.0E-7);
}

void Equilibrator::AddCondition(const string& type)
{
	// The new equilibrium condition
	EquilibriumCondition eCondition;
	
	if(type == "ChargeBalanced")
		eCondition = eSolver.ChargeBalanced();
	
	else 
	{ cerr << "Error: The equilibrium condition " << type << " is not valid." << endl; exit(0); }
	
	// Add the new equilibrium condition
	eConditions.push_back(eCondition);
}

void Equilibrator::AddCondition(const string& type, double value)
{
	// The new equilibrium condition
	EquilibriumCondition eCondition;
	
	if(type == "pH")
	{
		eCondition = eSolver.FixedAcidity(value);
		
		n[system["H[+]"]] = std::pow(10.0, -value);
	}
	
	else 
	{ cerr << "Error: The equilibrium condition " << type << " is not valid." << endl; exit(0); }
	
	// Add the new equilibrium condition
	eConditions.push_back(eCondition);
}

void Equilibrator::AddCondition(const string& type, const string& species, double value)
{
	// The new equilibrium condition
	EquilibriumCondition eCondition;
	
	// The global index of the species
	const Index ispecies = system[species];
	
	if(type == "FixedMass")
	{
		eCondition = eSolver.FixedMass(species, value);
		
		n[ispecies] = value;
	}
	
	else if(type == "FixedTotalMass")
	{
		eCondition = eSolver.FixedTotalMass(species, value);
		
		n[ispecies] = (value > 0.0) ? value : 1.0E-7;
	}
	
	else if(type == "FixedActivity")
	{
		eCondition = eSolver.FixedActivity(species, value);
		
		n[ispecies] = value;
	}
	
	else if(type == "FixedPartialPressure")
	{
		eCondition = eSolver.FixedPartialPressure(species, value);
		
		n[ispecies] = value;
	}
	
	else
	{ cerr << "Error: The equilibrium condition " << type << " is not valid." << endl; exit(0); }
	
	// Add the new equilibrium condition
	eConditions.push_back(eCondition);
}

void Equilibrator::Reset()
{
	// The number of species in the system
	const unsigned N = system.GetMultiphase().GetSpecies().size();
	
	// Clear the state of the composition of the system
	n = VectorXd::Constant(N, 1.0E-7);
	
	// Clear the equilibrium conditions provided so far
	eConditions.clear();
}

const SystemState Equilibrator::Equilibrate(double T, double P)
{
	eSolver.SetEquilibriumConditions(eConditions);
	
	eSolver.Solve(T, P, n);
	
	return SystemState(system, T, P, n);
}

